#include "AnaDef.h"
#include "TemplateFit.h"

//int main(int argc, char *argv[])
int templateFit(const Double_t kscalewid=1)
{
  //gErrorIgnoreLevel = 1001;

//
//  for(int ii=0; ii<argc; ii++){
//    printf("%d: %s\n", ii, argv[ii]);
//  }
//
  TRandom3 rndm(123);
  TH1D *hdata = new TH1D("hdata","",150,0,150);
  hdata->Sumw2();

  const Int_t ntmpl = 6;
  TH1D * htmpls[ntmpl];
  Double_t trueFrac[ntmpl];
  Double_t ntot=0;
  for(Int_t it=0; it<ntmpl; it++){
    htmpls[it] = new TH1D(Form("htmpl%d",it),"",150,0,150);
    htmpls[it]->Sumw2();

    const Int_t nsample = 1000*(it+1);
    trueFrac[it]=nsample;
    ntot+=nsample;

    const Double_t mean = 20+(it+1)*10;
    const Double_t wid = (0.5+(it+1)*0.2)*kscalewid;
    printf("it %d nsample %d mean %f wid %f\n", it, nsample, mean, wid);

    for(Int_t ip=0; ip<nsample; ip++){
      const Double_t var = rndm.Gaus(mean, wid);
      htmpls[it]->Fill(var);

      if(rndm.Rndm()<0.5){
        hdata->Fill(var);
      }
    }
    
    htmpls[it]->Scale(1./(htmpls[it]->Integral(0,10000)));
  } 
  for(Int_t it=0; it<ntmpl; it++){
    trueFrac[it] /= ntot;
  }

  hdata->Scale(1./(hdata->Integral(0,10000)));

  Double_t pars[ntmpl+2];
  for(Int_t ii=0; ii<ntmpl; ii++){
    pars[ii] = 1./ntmpl;
  }

  TemplateFit::Fit(hdata, ntmpl, htmpls, "draw", kFALSE, pars);
  TemplateFit::Fit(hdata, ntmpl, htmpls, "fit", 1, pars);

  printf("\nfit results: chi2 %.0f ndof %.0f \n", pars[ntmpl], pars[ntmpl+1]);
  for(Int_t ii=0; ii<ntmpl; ii++){
    printf("p%d: true %e fit %e fit/true %.2f\n", ii, trueFrac[ii], pars[ii], pars[ii]/trueFrac[ii]);
  }

  return 0;
}
